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To study the coexistence of two liquid states of water within one simulation box, we implement an 
equilibrium sedimentation method — which involves applying a gravitational field to the system and 
measuring/calculating the resulting density profile in equilibrium. We simulate a system of particles 
interacting via the ST2 potential, a model for water. We detect the coexistence of two liquid phases 
at low temperature. 
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I. INTRODUCTION 

The physics of the liquid state has been the subject 
of intense research activity. Novel approaches and novel 
apparatus have made possible the study of liquids under 
extreme thermodynamic conditions and in wide windows 
of space and time. Interesting and unexpected new phe- 
nomena have emerged as a result of a combined effort 
involving experiments, theory and simulations jjl ||. One 
of these is the possibility of a liquid-liquid (LL) transi- 
tion in one-component systems — in addition to the usual 
liquid-gas transition. Several liquids ||, |J, ^ and several 
models §, g, |, g, [k], [fl], [TJ, have been studied in 
detail, and it appears that the class of materials where 
a LL transition can be observed is larger than the class 
of tetrahedral liquids that were originally considered Jl4|] 
possible candidates for a LL transition. 

Water is one of the liquids that might possess a LL 
transition. Indeed, the first conjecture of a LL transition 
was based on a numerical study of the ST2 potential [^5) , 
a model designed to mimic the behavior of liquid water. 
In the case of ST2, the two coexisting phases differ in 
their local structure. The low density phase is formed by 
an open tetrahedrally-coordinated network of hydrogen 
bonds, while the high density phase has a more distorted 
network of hydrogen bonds. Recent theoretical work has 
shown that the interplay between local energy, entropy, 
and volume which may generate a LL transition can in 
principle be realized by spherically-symmetric potentials 
10. 

The evaluation of the P(p, T) equation of state (EOS) 
is key to test for phase coexistence (gas-liquid, liquid- 
crystal and LL), where P, p, and T denote the pressure, 
density, and temperature. The numerical calculation of 
P(p, To) requires the study of the model for a variety 
of state points. The coexistence between two phases, in 
the appropriate temperature window, appears as a region 
of density values where P is constant. In small-size nu- 



merical simulations is sometimes hard to observe phase 
transitions directly in one simulation box, in part because 
the free energy associated with creating an interface often 
stabilizes metastable phases pq. In these cases P(p,T) 
does not show any flat regionfll7|| . 

A different approach for studying in one single numer- 
ical simulation an entire isotherm has been proposed in 
Rcf. (and later exploited in the experimental study 
of colloidal systems and in the study of crystallization 
profiles Ref . [2(], |TJ ) . This approach simulates a 
semi-infinite tube in the presence of a very strong grav- 
itational field and measures the density profile in equi- 
librium. A simple inversion of the density profile allows 
the model EOS to be constructed. This idea has also 
been applied in the experimental study of the EOS of 
colloidal particles, by inversion of the measured sedimen- 
tation equilibrium profile. 

Here we apply the sedimentation profile method to 
ST2 water. We find that an interface separating two 
liquid states appears at low temperature, corresponding 
to the coexistence of two metastable liquid states of wa- 
ter within one simulation box, providing evidence for the 
presence of an LL transition. 



II. THEORY AND SIMULATION DETAILS 

We study N = 7680 rigid molecules of mass M in- 
teracting through the ST2 water potential, a rigid, non- 
polarizable, 5-site potential |Q that is able to repro- 
duce qualitatively the thermodynamic anomalies of liq- 
uid water. "ST2 water" is characterized, on cooling, by 
isobaric density maxima, increasing compressibility, in- 
creasing constant P specific heat, and evidence for a LL 
transition in the deeply supercooled regime (which is dif- 
ficult to probe experimentally due to spontaneous crys- 
tallization Jl5|). 

To implement the sedimentation profile method, we 
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FIG. 1: Schematic of the columnar box in a gravitational field 
g s . We partition the box into equally sized bins with height 
Az, where each bin represents a state point. The number 
of molecules in each bin is used to calculate the density and 
pressure at the center of the bin. 



use a column-shaped simulation box, semi-infinite along 
the z axis and with periodic boundary conditions along 
the x and y axes. The top of the box is left open, while 
the bottom is assumed to be a repulsive soft-sphere sur- 
face, generating a short-range force proportional to z -13 
acting on the molecule center of mass (Fig. [I]). The box 
width in the x-axis and y-axis directions is 3 nm, cor- 



responding to a bottom surface area of S x 



9 nm 2 . 



A strong gravitational field g s = 2 x I0 12 g is applied 
downward, in the z-axis direction, where g = 9.8 kg/m 2 
is the Earth's gravitational field. The value of the field 
controls the range of P-values accessed in the simula- 
tions. The pressure at the bottom of the column is 
P = g s NM/S xy ~ 500MPa. 

We perform the simulations using a multi-processor 
code on SGI Origin 2000, IBM SP, and IBM Regatta 
supercomputers. We choose a 1 fs time step and study 
different temperatures from T = 300 K down to T = 
230 K. For T=230 K, we simulate four different systems 
to better estimate error and reproducibility of the results. 
A long (~ 10ns) equilibration time precedes the actual 
calculation of the equilibrium density profile. To analyze 
equilibration we monitored the running average of the z- 
position of the center of mass. Production runs lasted at 
least 200 ps for the higher temperatures up to several ns 
per box when T = 230 K. 

We define a mass density field p(z) by averaging the 
density over bins with height Az = 1 nm to calculate the 
pressure field P(z) from 



P(z) 



p{z')dz'. 
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A parametric plot of P(z) vs. p(z) provides the EOS. 



FIG. 2: Density profiles for "ST2 water" for three different 
temperatures. At T=230 K a discontinuity in the density 
appears around 10 nm. In all boxes, the liquid-gas transition 
is at the state point at the maximum height shown. 



III. RESULTS AND DISCUSSION 

Figure ^| shows the equilibrium density profile along 
the z-axis for three different temperatures. Each symbol 
represents one bin, where the height of each bin Az = 
1 nm. For T — 300 K we see only one break in the density 
profile: the topmost points in the plots correspond to the 
gas-liquid interface. Above these points is a much less 
dense gas not shown on this density scale. For T = 250 K, 
there appears an inflection. For T = 230 K, we see a 
clear break in the density, associated with the interface 
between two different liquid states. 

Figure || shows the corresponding P{p) relations at the 
same three temperatures. We also show P{p) as evalu- 
ated previously using standard MD for cubic boxes with 
periodic boundary conditions for systems with N = 6 3 = 
216 (Ref. H) and N = 12 3 = 1768 (Ref.||). In 
Refs |24| and [p5| P was calculated with the standard 
virial relation Pq T 

Consistent with the data shown in Fig. we see a 
region of coexistence between two different phases at 
T = 230K, in agreement with the estimated location 
of the LL transition in Ref. ps| . Note that when using 
the sedimentation equilibration method, we see no un- 
physical loops in the equation of state, unlike cubic box 
simulations where the boundary conditions may artifi- 
cially stabilize metastablc states. 

To confirm that the two coexisting phases are both liq- 
uids, we calculate the mean square displacement (MSD) 
for different height values, both below and above the in- 
terface. We follow the evolution of each molecule for an 
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FIG. 3: EOS for the ST2 potential using the sedimenta- 
tion method shown as open circles. Open diamonds are the 
P(p,T) previously reported in Rehfcj], and the crosses are 
values previously reported in Ref. [|25[7\Ve corrected the data 
of Ref. J24] and Ref. |2f| to subtract the correction to P arising 
from the integration of the Lennard- Jones potential beyond 
the cut-off implemented in the simulation. Note the unphys- 
ical "loop" in the EOS at T = 230K which can be observed 
in the standard cubic-box simulations. 



average mean squared displacement smaller that Az, so 
that each height value can be unambiguously assigned 
to an average density value. Figure || shows that both 
phases are sufficiently diffusive. The low density phase 
has smaller diffusivity, in agreement with previous simu- 
lations of the density dependence of the dynamics p7J . 
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FIG. 4: (a) Mean square displacement of molecules for differ- 
ent density values within the column-shaped simulation box. 
Note that for values of the density both lower and higher than 
the coexisting phases, molecules diffuse with a finite diffusion 
constant, (b) Dependence of diffusion constant on density p 
across the interface. 



LL coexistence in the ST2 potential phase diagram is con- 
sistent with estimate of the "critical point" being located 
around T = 235 K, P = 250 MPa, and p = 1.05 g/cm 3 . 
Note that for most simple water models, the tempera- 
ture and pressure scales are shifted relative to real val- 
ues, thus they must be shifted to place them within an 
experimental context [^2[ . The critical point estimated 
experimentally is around T c — 230 K and P c = 50 MPa 
with p = 1.05 g/cm 3 pi]. 



IV. CONCLUSIONS 

We studied the EOS of "ST2 water" with the sedi- 
mentation profile method, and presented evidence that 
the EOS at T = 230-ftT shows a clear phase coexistence, 
between two phases which are both liquids. The observed 
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